library(tidyverse)
library(ggplot2)
library(ggpubr)
library(cowplot)
library(vegan)

Read in and Prep-process data

Read in Taxa Tables

find samples with 100% unknown and remove them

samples.unknown <- species_table %>% 
  melt () %>% 
  group_by(variable) %>% 
  summarise(sum=sum(value)) %>% 
  filter(sum == 0) %>% 
  pluck(1) %>% 
  as.vector()
No id variables; using all as measure variables
# Remove samples with no known taxa
species_table <- species_table[, -which(names(species_table) %in% samples.unknown)] 

Read in mapping table

remove neg controls; set sample.type factor

species_table <- species_table[,rownames(map_table)] # reorder to match metadata
species_table <- species_table[,-c(grep("Control", map_table$sample.type))] # remove control samples

genus_table <- genus_table[,rownames(map_table)]
genus_table <- genus_table[,-c(grep("Control", map_table$sample.type))]

map_table <- map_table[-c(grep("Control", map_table$sample.type)), ]
map_table$sample.type = factor(map_table$sample.type)

all(row.names(map_table) %in% colnames(species_table))
[1] TRUE

optional, remove follow up BAL

# species_table <- species_table[,rownames(map_table)] # rereorder to match metadata

# species_table <- species_table[,-c(grep("FALSE", map_table$is.baseline))]

#map_table <- map_table[-c(grep("FALSE", map_table$is.baseline)), ]

# set sample type as factor

#all(row.names(map_table) %in% colnames(species_table)) # Sanity check

Alpha diversity by sample type

compute shannon and simpson diversity metrics

diversity_vec = matrix(nrow = dim(species_table)[2], ncol = 2)
diversity_vec = as.data.frame(diversity_vec)
for (a in 1:dim(species_table)[2]) {
  diversity_vec[a,1] = diversity(species_table[,a], index = "shannon")
  diversity_vec[a,2] = diversity(species_table[,a], index = "simpson")
}
colnames(diversity_vec) = c("Shannon", "Simpson")

# add sample.type factor
diversity_vec$sample.type = map_table$sample.type
diversity_vec$sample.type = factor(diversity_vec$sample.type) # Optional: Add Levels

boxplots

ggplot(diversity_vec, aes(x = sample.type, y = Shannon, fill = sample.type)) +
  geom_boxplot(outlier.shape = NA, alpha = 0.7) +
  #geom_violin(alpha = 0.8) +
  geom_jitter(size = 1, width = 0.1, alpha = 0.35) +
  #scale_fill_manual(values=c("blue", "lightgreen", "yellow", "brown", "violet", "gray"),
  #                  labels = c("Toothbrush", "HMP-Oral", "HMP-Skin", "HMP-Gut", "HMP-Vaginal", "Building dust")) +
  theme_bw() +
  theme(legend.position = "none",
        axis.title.x = element_blank(),
        #axis.text.x = element_blank(), 
        axis.ticks.x = element_blank(), 
        axis.title.y = element_text(size = 14),
        axis.text.y = element_text(size = 14))

ANOVA stats

summary(aov(Shannon ~ sample.type, diversity_vec))
             Df Sum Sq Mean Sq F value Pr(>F)
sample.type   2   0.96  0.4790   0.896  0.409
Residuals   223 119.14  0.5343               
TukeyHSD(aov(Shannon ~ sample.type, diversity_vec))
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = Shannon ~ sample.type, data = diversity_vec)

$sample.type
                              diff        lwr       upr     p adj
Pneumonia-NonPneumonia -0.01538685 -0.4412036 0.4104299 0.9960005
Unknown-NonPneumonia    0.20087248 -0.3418635 0.7436084 0.6577194
Unknown-Pneumonia       0.21625932 -0.1650512 0.5975698 0.3755754

Beta diversity between sample types

# beta-diversity measure
beta <- vegdist(t(species_table), 'bray', binary = T)
beta.mat <- as.matrix(beta)

# projection
pcoa <- cmdscale(beta, k = 4, eig = TRUE)

# cleanup
ord <- as.data.frame(pcoa$points)
names(ord) <- c('pcoa1', 'pcoa2', 'pcoa3', 'pcoa4') ## rename coordinates

# add metadata
ord$Category = map_table$sample.type
ord$Outcome = map_table$binned_outcome
ord$Baseline = map_table$is.baseline
ord$Diagnosis = map_table$diagnosis.subtype

# Percent explained variation
eig <- eigenvals(pcoa)
eig.percent <- 100*head(eig/sum(eig))
eig.percent
[1] 27.568102 10.935027  9.314035  8.172572  5.591935  4.743136
eig.1 <- paste("PCo1 (", as.character(round(eig.percent[1], digits = 1)), "%)", sep = "")
eig.2 <-paste("PCo2 (", as.character(round(eig.percent[2], digits = 1)), "%)", sep = "")
## plot PCoA (FIGURE 1A)
ggplot(data = ord, aes(x = pcoa1, y = pcoa2, fill = Diagnosis)) +
  geom_point(size = 3, stroke = 1, shape=21) +
  theme_bw() +
  xlab(eig.1) +
  ylab(eig.2) +
  theme_q2r() +
  theme(axis.text = element_text(size=15, color = "black"),
        axis.title = element_text(size=16, color = "black"),
        legend.text = element_text(size=15, color = "black"),
        legend.title = element_text(size=15, color = "black")) +
  scale_fill_brewer(guide="legend", palette="Set1") +
  #scale_shape_manual(values=c(21))+#, 22, 23, 24, 25, 26)) +
  guides(fill=guide_legend(override.aes=list(shape=21))) +
   theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

ggsave("../../results/notebook_out/02A_MetaphlanPCoA.pdf", units="in", width=7, height=4.5)
Error in grDevices::pdf(file = filename, ..., version = version) : 
  cannot open file '../../results/notebook_out/02A_MetaphlanPCoA.pdf'

PERMANOVA Stats

# effect of sample type
beta <- vegdist(t(species_table), 'bray', binary = T)
adonis_out <- adonis2(beta ~ sample.type, data = map_table, permutations = 999)
adonis_out
Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = beta ~ sample.type, data = map_table, permutations = 999)
             Df SumOfSqs      R2      F Pr(>F)  
sample.type   2    1.364 0.01712 1.9416  0.017 *
Residual    223   78.308 0.98288                
Total       225   79.672 1.00000                
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Distance between samples

sample.type.1 <- map_table$sample.type[1]
sample.type.1
[1] Pneumonia
Levels: NonPneumonia Pneumonia Unknown
mean(as.matrix(vegdist(t(species_table[,grep(sample.type.1,map_table$sample.type)]), 
                       'jaccard', binary = T)))
[1] 0.8587207
sd(as.matrix(vegdist(t(species_table[,grep(sample.type.1,map_table$sample.type)]), 
                     'jaccard', binary = T)))/sqrt(34)
[1] 0.03923536

Heat Map Most Abundant Taxa

Get most abundant, prevalent species in ranked order

min_mean_proportion <- .00001
min_prevalence <- 20

species.ranked <- species_table %>% 
  rownames_to_column() %>% 
  melt() %>% 
  filter(value > 0) %>% # filter 0 values to get prevalence
  group_by(rowname) %>%
  summarise(mean = mean(value), n = n(), sum=sum(value)) %>% 
  arrange(desc(mean)) %>% 
  #filter(mean > min_mean_proportion) %>% # Filter by mean percent
  filter(n > min_prevalence) %>% # Filter by prevalence
  pluck(1)
Using rowname as id variables
length(species.ranked)
[1] 19

Order by species rank object and remove extra taxnomic info

species_table.heat <- species_table[species.ranked,] %>% 
  rownames_to_column() %>%
  separate(rowname, into = c("ExtraTaxa", "Taxa"), sep="s__") %>%
  select(!ExtraTaxa) %>% 
  column_to_rownames("Taxa") %>% 
  as.matrix
library(pheatmap)
annot.col <- map_table %>% select(sample.type, binned_outcome, is.baseline) 

paletteLength <- 100

myColors <- rev(colorRampPalette(rev(brewer.pal(n = 9, name ="Reds")))(paletteLength))
myColors <- rev(colorRampPalette(rev(brewer.pal(n = 9, name ="RdBu")))(paletteLength))

#pdf(file="../../results/notebook_out/02A_MetaphlanHeatPlot.pdf", width = 25, height = 20)
pheatmap(species_table.heat, 
         color = myColors,
         annotation_col = annot.col,
         angle_col = "45",
         show_colnames=FALSE,
         show_rownames = FALSE, scale="row")

#dev.off()

Genus

min_mean_proportion <- .00001
min_prevalence <- 20

genera.ranked <- genus_table %>% 
  rownames_to_column() %>% 
  melt() %>% 
  filter(value > 0) %>% # filter 0 values to get prevalence
  group_by(rowname) %>%
  summarise(mean = mean(value), n = n(), sum=sum(value)) %>% 
  arrange(desc(mean)) %>% 
  #filter(mean > min_mean_proportion) %>% # Filter by mean percent
  filter(n > min_prevalence) %>% # Filter by prevalence
  pluck(1)
Using rowname as id variables
length(genera.ranked)
[1] 13
# Order by species rank object and remove extra taxnomic info
genus_table.heat <- genus_table[genera.ranked,] %>% 
  rownames_to_column() %>%
  separate(rowname, into = c("ExtraTaxa", "Taxa"), sep="g__") %>%
  select(!ExtraTaxa) %>% 
  column_to_rownames("Taxa") %>% 
  as.matrix

Testing as phyloseq + Playground

LS0tCnRpdGxlOiAiMDJBX01ldGFQaGxhblRheG9ub21pY0FuYWx5c2lzIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgpgYGB7cn0KbGlicmFyeSh0aWR5dmVyc2UpCmxpYnJhcnkoZ2dwbG90MikKbGlicmFyeShnZ3B1YnIpCmxpYnJhcnkoY293cGxvdCkKbGlicmFyeSh2ZWdhbikKCmBgYAoKIyBSZWFkIGluIGFuZCBQcmVwLXByb2Nlc3MgZGF0YQojIyBSZWFkIGluIFRheGEgVGFibGVzCmBgYHtyfQoKdGF4YV90YWJsZSA9IHJlYWRfdGFibGUoJy4uLy4uL3Jlc3VsdHMvbWV0YXBobGFuX2Jvd3RpZV9vdXQvbWVyZ2VkX21ldGFwaGxhbl9wcm9maWxlLnRzdicsIHNraXAgPSAxKSAlPiUKICAgICAgICAgICAgICAgICAgICAgICAgc2VsZWN0KCFOQ0JJX3RheF9pZCkKCiMgUmVtb3ZlIC5tZXRhcGhsYW5fcHJvZmlsZSBmcm9tIGVuZHMgb2YgZmlsZXMgCmNvbG5hbWVzKHRheGFfdGFibGUpIDwtIHNhcHBseShzdHJzcGxpdChuYW1lcyh0YXhhX3RhYmxlKSwgIi5tZXRhcGhsYW5fcHJvZmlsZSIpLCBgW1tgLCAxKQoKZ2VudXNfdGFibGUgPC0gdGF4YV90YWJsZSAlPiUgZmlsdGVyKGdyZXBsKCJnX18iLCBjbGFkZV9uYW1lKSwgIWdyZXBsKCJzX18iLCBjbGFkZV9uYW1lKSkKZ2VudXNfdGFibGUgPC0gZ2VudXNfdGFibGUgJT4lIGNvbHVtbl90b19yb3duYW1lcyh2YXIgPSBjKCJjbGFkZV9uYW1lIikpCgpzcGVjaWVzX3RhYmxlIDwtIHRheGFfdGFibGUgJT4lIGZpbHRlcihncmVwbCgic19fIiwgY2xhZGVfbmFtZSkpCnNwZWNpZXNfdGFibGUgPC0gc3BlY2llc190YWJsZSAlPiUgY29sdW1uX3RvX3Jvd25hbWVzKHZhciA9IGMoImNsYWRlX25hbWUiKSkKYGBgCiMgZmluZCBzYW1wbGVzIHdpdGggMTAwJSB1bmtub3duIGFuZCByZW1vdmUgdGhlbQpgYGB7cn0Kc2FtcGxlcy51bmtub3duIDwtIHNwZWNpZXNfdGFibGUgJT4lIAogIG1lbHQgKCkgJT4lIAogIGdyb3VwX2J5KHZhcmlhYmxlKSAlPiUgCiAgc3VtbWFyaXNlKHN1bT1zdW0odmFsdWUpKSAlPiUgCiAgZmlsdGVyKHN1bSA9PSAwKSAlPiUgCiAgcGx1Y2soMSkgJT4lIAogIGFzLnZlY3RvcigpCgojIFJlbW92ZSBzYW1wbGVzIHdpdGggbm8ga25vd24gdGF4YQpzcGVjaWVzX3RhYmxlIDwtIHNwZWNpZXNfdGFibGVbLCAtd2hpY2gobmFtZXMoc3BlY2llc190YWJsZSkgJWluJSBzYW1wbGVzLnVua25vd24pXSAKCmBgYAoKIyMgUmVhZCBpbiBtYXBwaW5nIHRhYmxlCmBgYHtyfQptYXBfdGFibGUgPC0gcmVhZF90YWJsZSgnLi4vLi4vY29uZmlnL21hcF90YWJsZS50c3YnKSAlPiUKICAgICAgICAgICAgICAgICAgICAgICAgI2ZpbHRlcighc2FtcGxlICVpbiUgc2FtcGxlcy51bmtub3duKSAlPiUgIyByZW1vdmUgc2FtcGxlcyB3aXRoIHVua293bgogICAgICAgICAgICAgICAgICAgICAgICBjb2x1bW5fdG9fcm93bmFtZXModmFyID0gYygic2FtcGxlIikpCm1hcF90YWJsZSA8LSBtYXBfdGFibGVbLXdoaWNoKHJvd25hbWVzKG1hcF90YWJsZSkgJWluJSBzYW1wbGVzLnVua25vd24pLF0gCgpgYGAKCgojIyByZW1vdmUgbmVnIGNvbnRyb2xzOyBzZXQgc2FtcGxlLnR5cGUgZmFjdG9yCmBgYHtyfQpzcGVjaWVzX3RhYmxlIDwtIHNwZWNpZXNfdGFibGVbLHJvd25hbWVzKG1hcF90YWJsZSldICMgcmVvcmRlciB0byBtYXRjaCBtZXRhZGF0YQpzcGVjaWVzX3RhYmxlIDwtIHNwZWNpZXNfdGFibGVbLC1jKGdyZXAoIkNvbnRyb2wiLCBtYXBfdGFibGUkc2FtcGxlLnR5cGUpKV0gIyByZW1vdmUgY29udHJvbCBzYW1wbGVzCgpnZW51c190YWJsZSA8LSBnZW51c190YWJsZVsscm93bmFtZXMobWFwX3RhYmxlKV0KZ2VudXNfdGFibGUgPC0gZ2VudXNfdGFibGVbLC1jKGdyZXAoIkNvbnRyb2wiLCBtYXBfdGFibGUkc2FtcGxlLnR5cGUpKV0KCm1hcF90YWJsZSA8LSBtYXBfdGFibGVbLWMoZ3JlcCgiQ29udHJvbCIsIG1hcF90YWJsZSRzYW1wbGUudHlwZSkpLCBdCm1hcF90YWJsZSRzYW1wbGUudHlwZSA9IGZhY3RvcihtYXBfdGFibGUkc2FtcGxlLnR5cGUpCgphbGwocm93Lm5hbWVzKG1hcF90YWJsZSkgJWluJSBjb2xuYW1lcyhzcGVjaWVzX3RhYmxlKSkKCmBgYAoKIyMgb3B0aW9uYWwsIHJlbW92ZSBmb2xsb3cgdXAgQkFMCmBgYHtyfQojIHNwZWNpZXNfdGFibGUgPC0gc3BlY2llc190YWJsZVsscm93bmFtZXMobWFwX3RhYmxlKV0gIyByZXJlb3JkZXIgdG8gbWF0Y2ggbWV0YWRhdGEKCiMgc3BlY2llc190YWJsZSA8LSBzcGVjaWVzX3RhYmxlWywtYyhncmVwKCJGQUxTRSIsIG1hcF90YWJsZSRpcy5iYXNlbGluZSkpXQoKI21hcF90YWJsZSA8LSBtYXBfdGFibGVbLWMoZ3JlcCgiRkFMU0UiLCBtYXBfdGFibGUkaXMuYmFzZWxpbmUpKSwgXQoKIyBzZXQgc2FtcGxlIHR5cGUgYXMgZmFjdG9yCgojYWxsKHJvdy5uYW1lcyhtYXBfdGFibGUpICVpbiUgY29sbmFtZXMoc3BlY2llc190YWJsZSkpICMgU2FuaXR5IGNoZWNrCgpgYGAKCgoKIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjCiMgQWxwaGEgZGl2ZXJzaXR5IGJ5IHNhbXBsZSB0eXBlIAojIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMKCiMjIGNvbXB1dGUgc2hhbm5vbiBhbmQgc2ltcHNvbiBkaXZlcnNpdHkgbWV0cmljcwpgYGB7cn0KZGl2ZXJzaXR5X3ZlYyA9IG1hdHJpeChucm93ID0gZGltKHNwZWNpZXNfdGFibGUpWzJdLCBuY29sID0gMikKZGl2ZXJzaXR5X3ZlYyA9IGFzLmRhdGEuZnJhbWUoZGl2ZXJzaXR5X3ZlYykKZm9yIChhIGluIDE6ZGltKHNwZWNpZXNfdGFibGUpWzJdKSB7CiAgZGl2ZXJzaXR5X3ZlY1thLDFdID0gZGl2ZXJzaXR5KHNwZWNpZXNfdGFibGVbLGFdLCBpbmRleCA9ICJzaGFubm9uIikKICBkaXZlcnNpdHlfdmVjW2EsMl0gPSBkaXZlcnNpdHkoc3BlY2llc190YWJsZVssYV0sIGluZGV4ID0gInNpbXBzb24iKQp9CmNvbG5hbWVzKGRpdmVyc2l0eV92ZWMpID0gYygiU2hhbm5vbiIsICJTaW1wc29uIikKCiMgYWRkIHNhbXBsZS50eXBlIGZhY3RvcgpkaXZlcnNpdHlfdmVjJHNhbXBsZS50eXBlID0gbWFwX3RhYmxlJHNhbXBsZS50eXBlCmRpdmVyc2l0eV92ZWMkc2FtcGxlLnR5cGUgPSBmYWN0b3IoZGl2ZXJzaXR5X3ZlYyRzYW1wbGUudHlwZSkgIyBPcHRpb25hbDogQWRkIExldmVscwpgYGAKCgojIyBib3hwbG90cyAKYGBge3J9CmdncGxvdChkaXZlcnNpdHlfdmVjLCBhZXMoeCA9IHNhbXBsZS50eXBlLCB5ID0gU2hhbm5vbiwgZmlsbCA9IHNhbXBsZS50eXBlKSkgKwogIGdlb21fYm94cGxvdChvdXRsaWVyLnNoYXBlID0gTkEsIGFscGhhID0gMC43KSArCiAgI2dlb21fdmlvbGluKGFscGhhID0gMC44KSArCiAgZ2VvbV9qaXR0ZXIoc2l6ZSA9IDEsIHdpZHRoID0gMC4xLCBhbHBoYSA9IDAuMzUpICsKICAjc2NhbGVfZmlsbF9tYW51YWwodmFsdWVzPWMoImJsdWUiLCAibGlnaHRncmVlbiIsICJ5ZWxsb3ciLCAiYnJvd24iLCAidmlvbGV0IiwgImdyYXkiKSwKICAjICAgICAgICAgICAgICAgICAgbGFiZWxzID0gYygiVG9vdGhicnVzaCIsICJITVAtT3JhbCIsICJITVAtU2tpbiIsICJITVAtR3V0IiwgIkhNUC1WYWdpbmFsIiwgIkJ1aWxkaW5nIGR1c3QiKSkgKwogIHRoZW1lX2J3KCkgKwogIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICJub25lIiwKICAgICAgICBheGlzLnRpdGxlLnggPSBlbGVtZW50X2JsYW5rKCksCiAgICAgICAgI2F4aXMudGV4dC54ID0gZWxlbWVudF9ibGFuaygpLCAKICAgICAgICBheGlzLnRpY2tzLnggPSBlbGVtZW50X2JsYW5rKCksIAogICAgICAgIGF4aXMudGl0bGUueSA9IGVsZW1lbnRfdGV4dChzaXplID0gMTQpLAogICAgICAgIGF4aXMudGV4dC55ID0gZWxlbWVudF90ZXh0KHNpemUgPSAxNCkpCmBgYAojIyBBTk9WQSBzdGF0cwpgYGB7cn0Kc3VtbWFyeShhb3YoU2hhbm5vbiB+IHNhbXBsZS50eXBlLCBkaXZlcnNpdHlfdmVjKSkKVHVrZXlIU0QoYW92KFNoYW5ub24gfiBzYW1wbGUudHlwZSwgZGl2ZXJzaXR5X3ZlYykpCmBgYAojIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMKIyBCZXRhIGRpdmVyc2l0eSBiZXR3ZWVuIHNhbXBsZSB0eXBlcyAKIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjCgpgYGB7cn0KIyBiZXRhLWRpdmVyc2l0eSBtZWFzdXJlCmJldGEgPC0gdmVnZGlzdCh0KHNwZWNpZXNfdGFibGUpLCAnYnJheScsIGJpbmFyeSA9IFQpCmJldGEubWF0IDwtIGFzLm1hdHJpeChiZXRhKQoKIyBwcm9qZWN0aW9uCnBjb2EgPC0gY21kc2NhbGUoYmV0YSwgayA9IDQsIGVpZyA9IFRSVUUpCgojIGNsZWFudXAKb3JkIDwtIGFzLmRhdGEuZnJhbWUocGNvYSRwb2ludHMpCm5hbWVzKG9yZCkgPC0gYygncGNvYTEnLCAncGNvYTInLCAncGNvYTMnLCAncGNvYTQnKSAjIyByZW5hbWUgY29vcmRpbmF0ZXMKCiMgYWRkIG1ldGFkYXRhCm9yZCRDYXRlZ29yeSA9IG1hcF90YWJsZSRzYW1wbGUudHlwZQpvcmQkT3V0Y29tZSA9IG1hcF90YWJsZSRiaW5uZWRfb3V0Y29tZQpvcmQkQmFzZWxpbmUgPSBtYXBfdGFibGUkaXMuYmFzZWxpbmUKb3JkJERpYWdub3NpcyA9IG1hcF90YWJsZSRkaWFnbm9zaXMuc3VidHlwZQoKIyBQZXJjZW50IGV4cGxhaW5lZCB2YXJpYXRpb24KZWlnIDwtIGVpZ2VudmFscyhwY29hKQplaWcucGVyY2VudCA8LSAxMDAqaGVhZChlaWcvc3VtKGVpZykpCmVpZy5wZXJjZW50CmVpZy4xIDwtIHBhc3RlKCJQQ28xICgiLCBhcy5jaGFyYWN0ZXIocm91bmQoZWlnLnBlcmNlbnRbMV0sIGRpZ2l0cyA9IDEpKSwgIiUpIiwgc2VwID0gIiIpCmVpZy4yIDwtcGFzdGUoIlBDbzIgKCIsIGFzLmNoYXJhY3Rlcihyb3VuZChlaWcucGVyY2VudFsyXSwgZGlnaXRzID0gMSkpLCAiJSkiLCBzZXAgPSAiIikKCmBgYAoKYGBge3J9CiMjIHBsb3QgUENvQSAoRklHVVJFIDFBKQpnZ3Bsb3QoZGF0YSA9IG9yZCwgYWVzKHggPSBwY29hMSwgeSA9IHBjb2EyLCBmaWxsID0gRGlhZ25vc2lzKSkgKwogIGdlb21fcG9pbnQoc2l6ZSA9IDMsIHN0cm9rZSA9IDEsIHNoYXBlPTIxKSArCiAgdGhlbWVfYncoKSArCiAgeGxhYihlaWcuMSkgKwogIHlsYWIoZWlnLjIpICsKICB0aGVtZV9xMnIoKSArCiAgdGhlbWUoYXhpcy50ZXh0ID0gZWxlbWVudF90ZXh0KHNpemU9MTUsIGNvbG9yID0gImJsYWNrIiksCiAgICAgICAgYXhpcy50aXRsZSA9IGVsZW1lbnRfdGV4dChzaXplPTE2LCBjb2xvciA9ICJibGFjayIpLAogICAgICAgIGxlZ2VuZC50ZXh0ID0gZWxlbWVudF90ZXh0KHNpemU9MTUsIGNvbG9yID0gImJsYWNrIiksCiAgICAgICAgbGVnZW5kLnRpdGxlID0gZWxlbWVudF90ZXh0KHNpemU9MTUsIGNvbG9yID0gImJsYWNrIikpICsKICBzY2FsZV9maWxsX2JyZXdlcihndWlkZT0ibGVnZW5kIiwgcGFsZXR0ZT0iU2V0MSIpICsKICAjc2NhbGVfc2hhcGVfbWFudWFsKHZhbHVlcz1jKDIxKSkrIywgMjIsIDIzLCAyNCwgMjUsIDI2KSkgKwogIGd1aWRlcyhmaWxsPWd1aWRlX2xlZ2VuZChvdmVycmlkZS5hZXM9bGlzdChzaGFwZT0yMSkpKSArCiAgIHRoZW1lKHBhbmVsLmdyaWQubWFqb3IgPSBlbGVtZW50X2JsYW5rKCksIHBhbmVsLmdyaWQubWlub3IgPSBlbGVtZW50X2JsYW5rKCkpIApnZ3NhdmUoIi4uLy4uL3Jlc3VsdHMvbm90ZWJvb2tfb3V0LzAyQV9NZXRhcGhsYW5QQ29BLnBkZiIsIHVuaXRzPSJpbiIsIHdpZHRoPTcsIGhlaWdodD00LjUpCmBgYAoKIyMgUEVSTUFOT1ZBIFN0YXRzCmBgYHtyfQojIGVmZmVjdCBvZiBzYW1wbGUgdHlwZQpiZXRhIDwtIHZlZ2Rpc3QodChzcGVjaWVzX3RhYmxlKSwgJ2JyYXknLCBiaW5hcnkgPSBUKQphZG9uaXNfb3V0IDwtIGFkb25pczIoYmV0YSB+IHNhbXBsZS50eXBlLCBkYXRhID0gbWFwX3RhYmxlLCBwZXJtdXRhdGlvbnMgPSA5OTkpCmFkb25pc19vdXQKYGBgCgojIyBEaXN0YW5jZSBiZXR3ZWVuIHNhbXBsZXMKYGBge3J9CnNhbXBsZS50eXBlLjEgPC0gbWFwX3RhYmxlJHNhbXBsZS50eXBlWzFdCnNhbXBsZS50eXBlLjEKCm1lYW4oYXMubWF0cml4KHZlZ2Rpc3QodChzcGVjaWVzX3RhYmxlWyxncmVwKHNhbXBsZS50eXBlLjEsbWFwX3RhYmxlJHNhbXBsZS50eXBlKV0pLCAKICAgICAgICAgICAgICAgICAgICAgICAnamFjY2FyZCcsIGJpbmFyeSA9IFQpKSkKc2QoYXMubWF0cml4KHZlZ2Rpc3QodChzcGVjaWVzX3RhYmxlWyxncmVwKHNhbXBsZS50eXBlLjEsbWFwX3RhYmxlJHNhbXBsZS50eXBlKV0pLCAKICAgICAgICAgICAgICAgICAgICAgJ2phY2NhcmQnLCBiaW5hcnkgPSBUKSkpL3NxcnQoMzQpCmBgYAoKIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjCiMgSGVhdCBNYXAgTW9zdCBBYnVuZGFudCBUYXhhCiMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIwoKIyMgR2V0IG1vc3QgYWJ1bmRhbnQsIHByZXZhbGVudCBzcGVjaWVzIGluIHJhbmtlZCBvcmRlcgpgYGB7cn0KbWluX21lYW5fcHJvcG9ydGlvbiA8LSAuMDAwMDEKbWluX3ByZXZhbGVuY2UgPC0gMjAKCnNwZWNpZXMucmFua2VkIDwtIHNwZWNpZXNfdGFibGUgJT4lIAogIHJvd25hbWVzX3RvX2NvbHVtbigpICU+JSAKICBtZWx0KCkgJT4lIAogIGZpbHRlcih2YWx1ZSA+IDApICU+JSAjIGZpbHRlciAwIHZhbHVlcyB0byBnZXQgcHJldmFsZW5jZQogIGdyb3VwX2J5KHJvd25hbWUpICU+JQogIHN1bW1hcmlzZShtZWFuID0gbWVhbih2YWx1ZSksIG4gPSBuKCksIHN1bT1zdW0odmFsdWUpKSAlPiUgCiAgYXJyYW5nZShkZXNjKG1lYW4pKSAlPiUgCiAgI2ZpbHRlcihtZWFuID4gbWluX21lYW5fcHJvcG9ydGlvbikgJT4lICMgRmlsdGVyIGJ5IG1lYW4gcGVyY2VudAogIGZpbHRlcihuID4gbWluX3ByZXZhbGVuY2UpICU+JSAjIEZpbHRlciBieSBwcmV2YWxlbmNlCiAgcGx1Y2soMSkKCmxlbmd0aChzcGVjaWVzLnJhbmtlZCkKYGBgCgojIyBPcmRlciBieSBzcGVjaWVzIHJhbmsgb2JqZWN0IGFuZCByZW1vdmUgZXh0cmEgdGF4bm9taWMgaW5mbwpgYGB7cn0Kc3BlY2llc190YWJsZS5oZWF0IDwtIHNwZWNpZXNfdGFibGVbc3BlY2llcy5yYW5rZWQsXSAlPiUgCiAgcm93bmFtZXNfdG9fY29sdW1uKCkgJT4lCiAgc2VwYXJhdGUocm93bmFtZSwgaW50byA9IGMoIkV4dHJhVGF4YSIsICJUYXhhIiksIHNlcD0ic19fIikgJT4lCiAgc2VsZWN0KCFFeHRyYVRheGEpICU+JSAKICBjb2x1bW5fdG9fcm93bmFtZXMoIlRheGEiKSAlPiUgCiAgYXMubWF0cml4CmBgYAoKYGBge3J9CmxpYnJhcnkocGhlYXRtYXApCmFubm90LmNvbCA8LSBtYXBfdGFibGUgJT4lIHNlbGVjdChzYW1wbGUudHlwZSwgYmlubmVkX291dGNvbWUsIGlzLmJhc2VsaW5lKSAKCnBhbGV0dGVMZW5ndGggPC0gMTAwCgpteUNvbG9ycyA8LSByZXYoY29sb3JSYW1wUGFsZXR0ZShyZXYoYnJld2VyLnBhbChuID0gOSwgbmFtZSA9IlJlZHMiKSkpKHBhbGV0dGVMZW5ndGgpKQpteUNvbG9ycyA8LSByZXYoY29sb3JSYW1wUGFsZXR0ZShyZXYoYnJld2VyLnBhbChuID0gOSwgbmFtZSA9IlJkQnUiKSkpKHBhbGV0dGVMZW5ndGgpKQoKI3BkZihmaWxlPSIuLi8uLi9yZXN1bHRzL25vdGVib29rX291dC8wMkFfTWV0YXBobGFuSGVhdFBsb3QucGRmIiwgd2lkdGggPSAyNSwgaGVpZ2h0ID0gMjApCnBoZWF0bWFwKHNwZWNpZXNfdGFibGUuaGVhdCwgCiAgICAgICAgIGNvbG9yID0gbXlDb2xvcnMsCiAgICAgICAgIGFubm90YXRpb25fY29sID0gYW5ub3QuY29sLAogICAgICAgICBhbmdsZV9jb2wgPSAiNDUiLAogICAgICAgICBzaG93X2NvbG5hbWVzPUZBTFNFLAogICAgICAgICBzaG93X3Jvd25hbWVzID0gRkFMU0UsIHNjYWxlPSJyb3ciKQojZGV2Lm9mZigpCmBgYAoKCiMjIEdlbnVzCmBgYHtyfQptaW5fbWVhbl9wcm9wb3J0aW9uIDwtIC4wMDAwMQptaW5fcHJldmFsZW5jZSA8LSAyMAoKZ2VuZXJhLnJhbmtlZCA8LSBnZW51c190YWJsZSAlPiUgCiAgcm93bmFtZXNfdG9fY29sdW1uKCkgJT4lIAogIG1lbHQoKSAlPiUgCiAgZmlsdGVyKHZhbHVlID4gMCkgJT4lICMgZmlsdGVyIDAgdmFsdWVzIHRvIGdldCBwcmV2YWxlbmNlCiAgZ3JvdXBfYnkocm93bmFtZSkgJT4lCiAgc3VtbWFyaXNlKG1lYW4gPSBtZWFuKHZhbHVlKSwgbiA9IG4oKSwgc3VtPXN1bSh2YWx1ZSkpICU+JSAKICBhcnJhbmdlKGRlc2MobWVhbikpICU+JSAKICAjZmlsdGVyKG1lYW4gPiBtaW5fbWVhbl9wcm9wb3J0aW9uKSAlPiUgIyBGaWx0ZXIgYnkgbWVhbiBwZXJjZW50CiAgZmlsdGVyKG4gPiBtaW5fcHJldmFsZW5jZSkgJT4lICMgRmlsdGVyIGJ5IHByZXZhbGVuY2UKICBwbHVjaygxKQoKbGVuZ3RoKGdlbmVyYS5yYW5rZWQpCmBgYApgYGB7cn0KIyBPcmRlciBieSBzcGVjaWVzIHJhbmsgb2JqZWN0IGFuZCByZW1vdmUgZXh0cmEgdGF4bm9taWMgaW5mbwpnZW51c190YWJsZS5oZWF0IDwtIGdlbnVzX3RhYmxlW2dlbmVyYS5yYW5rZWQsXSAlPiUgCiAgcm93bmFtZXNfdG9fY29sdW1uKCkgJT4lCiAgc2VwYXJhdGUocm93bmFtZSwgaW50byA9IGMoIkV4dHJhVGF4YSIsICJUYXhhIiksIHNlcD0iZ19fIikgJT4lCiAgc2VsZWN0KCFFeHRyYVRheGEpICU+JSAKICBjb2x1bW5fdG9fcm93bmFtZXMoIlRheGEiKSAlPiUgCiAgYXMubWF0cml4CmBgYAoKYGBge3J9CmxpYnJhcnkocGhlYXRtYXApCmFubm90LmNvbCA8LSBtYXBfdGFibGUgJT4lIHNlbGVjdChzYW1wbGUudHlwZSwgYmlubmVkX291dGNvbWUsIGlzLmJhc2VsaW5lKSAKCnBhbGV0dGVMZW5ndGggPC0gNTAKCm15Q29sb3JzIDwtIHJldihjb2xvclJhbXBQYWxldHRlKHJldihicmV3ZXIucGFsKG4gPSA5LCBuYW1lID0iUmVkcyIpKSkocGFsZXR0ZUxlbmd0aCkpCm15Q29sb3JzIDwtIHJldihjb2xvclJhbXBQYWxldHRlKHJldihicmV3ZXIucGFsKG4gPSA5LCBuYW1lID0iUmRCdSIpKSkocGFsZXR0ZUxlbmd0aCkpCgojcGRmKGZpbGU9Ii4uLy4uL3Jlc3VsdHMvbm90ZWJvb2tfb3V0LzAyQV9NZXRhcGhsYW5IZWF0UGxvdEdlbmVyYS5wZGYiLCB3aWR0aCA9IDI1LCBoZWlnaHQgPSAyMCkKcGhlYXRtYXAoZ2VudXNfdGFibGUuaGVhdCwgCiAgICAgICAgIGNvbG9yID0gbXlDb2xvcnMsCiAgICAgICAgIGFubm90YXRpb25fY29sID0gYW5ub3QuY29sLAogICAgICAgICBhbmdsZV9jb2wgPSAiNDUiLAogICAgICAgICBzaG93X2NvbG5hbWVzPUZBTFNFLAogICAgICAgICBzaG93X3Jvd25hbWVzID0gRkFMU0UsIHNjYWxlPSJyb3ciKQojZGV2Lm9mZigpCmBgYAoKCiMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIwojIFRlc3RpbmcgYXMgcGh5bG9zZXEgKyBQbGF5Z3JvdW5kCiMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIwoKYGBge3J9CmxpYnJhcnkodGlkeXZlcnNlKTsgcGFja2FnZVZlcnNpb24oInRpZHl2ZXJzZSIpICAgICAjdmVyc2lvbjoxLjMuMCAKbGlicmFyeShwaHlsb3NlcSk7IHBhY2thZ2VWZXJzaW9uKCJwaHlsb3NlcSIpICAgICAgICN2ZXJzaW9uOjEuMzIuMAoKCiNCdWlsZGluZyBNZXRhUGhsQW4gc3BlY2llcyBhYnVuZGFuY2UgcHMgb2JqZWN0CnNfYWJ1bmQgPC0gcmVhZF90c3YoIi4uLy4uL3Jlc3VsdHMvbWV0YXBobGFuX2Jvd3RpZV9vdXQvbWVyZ2VkX21ldGFwaGxhbl9wcm9maWxlX2FsbC50c3YiKSAlPiUgc2VsZWN0KCFOQ0JJX3RheF9pZCkKCnNfdGF4X3RhYiA8LSBzX2FidW5kICU+JQogIGRwbHlyOjpyZW5hbWUoInRheG9ub215IiA9ICJjbGFkZV9uYW1lIikgJT4lCiAgZHBseXI6OnNlbGVjdCh0YXhvbm9teSkgJT4lCiAgdGlkeXI6OnNlcGFyYXRlKHRheG9ub215LCBpbnRvID0gYygiS2luZ2RvbSIsICJQaHlsdW0iLCAiQ2xhc3MiLCAiT3JkZXIiLCAiRmFtaWx5IiwgIkdlbnVzIiwgIlNwZWNpZXMiKSwgc2VwID0gIlxcfCIpICU+JQogIGRwbHlyOjptdXRhdGUoc3BlY19yb3cgPSBTcGVjaWVzKSAlPiUKICB0aWJibGU6OmNvbHVtbl90b19yb3duYW1lcyh2YXIgPSAic3BlY19yb3ciKQoKc19vdHVfdGFiIDwtIHNfYWJ1bmQgJT4lCiAgZHBseXI6OnJlbmFtZSgidGF4b25vbXkiID0gImNsYWRlX25hbWUiKSAlPiUKICB0aWR5cjo6c2VwYXJhdGUodGF4b25vbXksIGludG8gPSBjKCJLaW5nZG9tIiwgIlBoeWx1bSIsICJDbGFzcyIsICJPcmRlciIsICJGYW1pbHkiLCAiR2VudXMiLCAiU3BlY2llcyIpLCBzZXAgPSAiXFx8IikgJT4lCiAgZHBseXI6OnNlbGVjdCgtS2luZ2RvbSwgLVBoeWx1bSwgLUNsYXNzLCAtT3JkZXIsIC1GYW1pbHksIC1HZW51cykgJT4lCiAgdGliYmxlOjpjb2x1bW5fdG9fcm93bmFtZXModmFyID0gIlNwZWNpZXMiKQoKbmFtZXMoc19vdHVfdGFiKSA8LSBnc3ViKG5hbWVzKHNfb3R1X3RhYiksIHBhdHRlcm4gPSAiLm1ldGFwaGxhbl9wcm9maWxlIiwgcmVwbGFjZW1lbnQgPSAiIikgCgpoZWFkKGNvbFN1bXMoc19vdHVfdGFiKSkKc19vdHVfdGFiIDwtIHNfb3R1X3RhYiAvIDEwMCAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICNjb252ZXJ0IHRvIHByb3BvcnRpb24gd2l0aCB1bml0IHN1bSBvZiAxCmhlYWQoY29sU3VtcyhzX290dV90YWIpKQoKc19tZXRhIDwtIGRhdGEuZnJhbWUoc2VxX2lkID0gbmFtZXMoc19vdHVfdGFiKSkKc19tZXRhIDwtIHNfbWV0YSAlPiUKICBkcGx5cjo6bXV0YXRlKHNhbXBsZU5hbWVzX3JvdyA9IHNlcV9pZCkgJT4lCiAgdGliYmxlOjpjb2x1bW5fdG9fcm93bmFtZXModmFyID0gInNhbXBsZU5hbWVzX3JvdyIpCgoocHNfbXBhM19zcGVjaWVzIDwtIHBoeWxvc2VxKHNhbXBsZV9kYXRhKHNfbWV0YSksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgb3R1X3RhYmxlKHNfb3R1X3RhYiwgdGF4YV9hcmVfcm93cyA9IFRSVUUpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIHRheF90YWJsZShhcy5tYXRyaXgoc190YXhfdGFiKSkpKQpsaWJyYXJ5KHZlZ2FuKQpkZWNvc3RhbmQoc19vdHVfdGFiLCBtZXRob2Q9Im5vcm1hbGl6ZSIsIE1BUkdJTj0xKQpgYGAK